function [spikes Gampa]=izhikevich_fs(stim_conductance_exc, stim_conductance_inh)

C=20;vr=-55;vt=-40;k=1;
a=0.15;b=8;c=-55;d=200;
vpeak=25;

T=1000;

dt_ms=1;


spikes = 0;
n=round(T/dt_ms);
v=vr*ones(1,n); u=0*v;
si= stim_conductance_inh .*[ones(1,n)];
se= stim_conductance_exc .*[ones(1,n)];

% model conductances.
tau_ampa  = exp(-dt_ms/5);
tau_nmda  = exp(-dt_ms/150);
tau_gabaa = exp(-dt_ms/6);
tau_gabab = exp(-dt_ms/150);

E_AMPA	= 0.0;
E_NMDA	= 0.0;
E_GABAa	= -70;
E_GABAb = -90;
Ggabaa = 0;Ggabab = 0; Gampa = 0;Gnmda = 0;

Ggabaa=zeros(1,n);
Ggabab=zeros(1,n);
Gampa=zeros(1,n);
Gnmda=zeros(1,n);
I=zeros(1,n);

gain_ampa=1;
gain_nmda=.5;
gain_gabaa=1;
gain_gabab=.1;
% end model conductances.


% stsp parameters:

tau_F1=exp(-dt_ms/1000);
tau_D1=exp(-dt_ms/800);
U=0.5;

stsp = U*ones(1,n);
stf= zeros(1,n);
stdep = zeros(1,n);

% end stsp paramters

for i = 1:n-1
    
    % conductances from synapses.
    Ggabaa(i+1) = Ggabaa(i)*tau_gabaa + dt_ms*stsp(i)*gain_gabaa*si(i);
    Ggabab(i+1) = Ggabab(i)*tau_gabab + dt_ms*stsp(i)*gain_gabab*si(i);
    
    Gampa(i+1) = Gampa(i)*tau_ampa + dt_ms*stsp(i)*gain_ampa*se(i);
    Gnmda(i+1) = Gnmda(i)*tau_nmda + dt_ms*stsp(i)*gain_nmda*se(i);

    xx = (-80-v(i))/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
    
     
    %g= ( Gampa(i+1)       +NMDAgate.*Gnmda(i+1)        + Ggabaa(i+1)         + Ggabab(i+1));
	%E= ( Gampa(i+1)*E_AMPA+NMDAgate.*Gnmda(i+1)*E_NMDA + Ggabaa(i+1)*E_GABAa + Ggabab(i+1)*E_GABAb);
    %I(i) = v(i).*g-E;
    I(i) = (v(i)-E_AMPA)*Gampa(i+1) + (v(i)-E_NMDA)*NMDAgate*Gnmda(i+1)+ ...
        (v(i)-E_GABAa)*Ggabaa(i+1)+(v(i)-E_GABAb)*Ggabab(i+1);
    
    
    v(i+1)=v(i)+.5*dt_ms*(k*(v(i)-vr)*(v(i)-vt)-u(i)-I(i))/C;
    v(i+1)=v(i+1)+.5*dt_ms*(k*(v(i+1)-vr)*(v(i+1)-vt)-u(i)-I(i))/C;
    u(i+1)=u(i)+dt_ms*a*(b*(v(i)-vr)-u(i));
    
    % Short term synaptic plasticity.
    stf(i+1) = stf(i)*tau_F1;
	stdep(i+1) = stdep(i)*tau_D1;
    stsp(i+1) = (1-stdep(i))*(U+stf(i)); % By default, keep the value of the last spike?
    % stsp
    if se(i) > 0.001
        % pretend like this means a presynaptic spike.
        stdep(i+1)=stdep(i)+stsp(i+1);
        stf(i+1)=stf(i+1)+U*(1-U-stf(i+1));
    end
    
    
    if v(i+1)>vpeak
        v(i)=vpeak;
        v(i+1)=c;
        u(i+1)=u(i+1)+d;
        spikes = spikes + 1;
        
        
    end
end
figure(1)
subplot(4,1,1)
plot(dt_ms*(1:n),[v' u']);
subplot(4,1,2)
plot(dt_ms*(1:n),[Gampa' Gnmda' Ggabaa' Ggabab']);
title('conductances');
subplot(4,1,3);
plot(dt_ms*(1:n),I);
title('Isyn (pA)');
subplot(4,1,4);
plot(dt_ms*(1:n),[stsp' stdep' stf' ]);
title('stsp');
legend('stsp','stdep','stf');
drawnow

